SUBROUTINE readtherm(isp)
	
	USE ckvar

	IMPLICIT NONE
	INTEGER :: isp
	REAL :: mwi

	INTEGER :: io,i,j
	INTEGER :: nelem(3)
	CHARACTER(8) :: tzz,txx
	CHARACTER(3) :: elem(3)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! READ thermo properties from therm.dat & open file !!
	open(unit=10,file='UCSD/therm.dat')
	
	io=0
	i=0
	j=0
	nelem=0
	DO
	i=i+1
	IF(io.ne.0.or.tzz.eq.'END') EXIT
	if(i.eq.1) then
	read(10,*,IOSTAT=io) tzz
	read(10,*) cktmin,cktmed,cktmax
	i=2
	endif
	IF(i.ge.2) THEN
	read(10,100,IOSTAT=io) tzz,txx,elem(1),nelem(1) &
				      ,elem(2),nelem(2) &
				      ,elem(3),nelem(3)
	read(10,'(5(1pe15.8))',IOSTAT=io) (thcef(isp,j),j=1,5)
	read(10,'(5(1pe15.8))',IOSTAT=io) (thcef(isp,j),j=6,10)
	read(10,'(4(1pe15.8))',IOSTAT=io) (thcef(isp,j),j=11,14)
	if(tzz.eq.spc(isp)) then
!	write(*,*) 'Find ',spc(isp),'!!'
	EXIT
	endif
	ENDIF
	ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Calculate molecular weight from mw.dat !!
	mwspc(isp)=0.0
	DO i=1,3
	IF(nelem(i).ne.0) then
	call readmw(elem(i),mwi)
	mwspc(isp)=mwspc(isp)+mwi*nelem(i)
	ENDIF
	ENDDO
!! End of calculating molecular weight    !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!	write(*,*) tzz
!	write(*,*) (elem(j),j=1,3)
!	write(*,*) (nelem(j),j=1,3)
	
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Calculate thermodyanamic properties w/ temp and cef(j) !!
!! thp(1)=Cp(T)/R					   !!
!! thp(2)=H^o(T)/RT					   !!
!! thp(3)=S^o(T)/R					   !!
!! Cp(T): heat capacity [kcal/mol-K]			   !!
!! H^o(T): standard state enthalpy [kcal/mol]		   !!
!! S^o(T): standard state entropy [kcal/mol-K]		   !!
!! R=0.001987 kcal/mol-K				   !! 
!	Ru=0.001987	!! Universal gas constant
!	if(temp.ge.tmin.and.temp.le.tmed) then
!	write(*,*) 'Low Temperature'
!	i=1	!! for low temperature
!	elseif(temp.gt.tmed.and.temp.le.tmax) then
!	write(*,*) 'High Temperature'
!	i=8
!	else
!	write(*,*) 'Temp must be within',tmin,'and',tmax
!	i=1
!	endif
!
!	thp(1)=thcef(i) + thcef(i+1)*temp + thcef(i+2)*temp**2 + &
!	       thcef(i+3)*temp**3 + cef(i+4)*temp**4
!	thp(2)=thcef(i) + thcef(i+1)*temp/2 + thcef(i+2)*temp**2/3 + &
!	       thcef(i+3)*temp**3/4 + thcef(i+4)*temp**4/5 + &
!	       thcef(i+5)/temp
!	thp(3)=thcef(i)*log(temp) + thcef(i+1)*temp + &
!	       thcef(i+2)*temp**2/2 + thcef(i+3)*temp**3/3 + &
!	       thcef(i+4)*temp**4/4 + thcef(i+6)
!
!	write(*,*) 'Temp=',temp
!	write(*,*) 'Cp(T)=',thp(1)*Ru
!	write(*,*) 'H^o(T)=',thp(2)*Ru*temp
!	write(*,*) 'S^o(T)=',thp(3)*Ru
!!							   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	close(10)
  100	format(A8,A16,3(A3,I2))
END SUBROUTINE readtherm
